#include <iostream>
#include "Matrix.h"
#include <cassert>
#include <cmath>

Matrix::Matrix(const Matrix& matrix)
{
	setNumCol(matrix.getNumCol());
	setNumRow(matrix.getNumRow());
	
	array = new double[matrix.getSize()];
	double *pm = matrix.getpArray();
	
	for(int i = 0; i < getSize(); i++)
		array[i] = pm[i];	
}


Matrix& Matrix::operator=(const Matrix& A)
{
	if(this == &A)
		return *this;
	
	int t = getSize();
	int a = A.getSize();
	
	if(t < a)
	{
		delete[] array;	
		array = new double[a];
		setNumCol(A.getNumCol());
		setNumRow(A.getNumRow());
	}
	
	double *pA = A.getpArray();
		
	for(int i = 0; i < getSize(); i++)
		array[i] = pA[i];

	
	return *this;
}

Matrix& Matrix::operator+=(const Matrix& B)
{
	double *pB = B.getpArray();
	for(int i = 0; i < getSize(); i++)
		array[i] += pB[i];
	
	return *this;
}

Matrix& Matrix::operator-=(const Matrix& B)
{
	double *pB = B.getpArray();
	for(int i = 0; i < getSize(); i++)
		array[i] -= pB[i];
	
	return *this;
}

Matrix& Matrix::operator*=(const Matrix& B)
{
	double c[getNumRow()*B.getNumCol()];
	double *b = B.getpArray();
	
	for(int i = 0, ii = 0, p=0; p < getNumRow(); i += getNumCol(), p++, ii += B.getNumCol())
		for(int j = 0; j < B.getNumCol(); j++)
		{
			int index = ((row == 1) || (B.getNumCol() == 1) ? p+j : ii+j);
			c[index] = 0;
			for(int k = 0, t = 0; k < B.getNumRow(); k++, t += B.getNumCol())
				c[index] += array[i+k]*b[t+j];
		}
				
		setNumCol(B.getNumCol());
		if(getSize() < getNumRow()*B.getNumCol())		
		if(array)
		{
			delete[] array;
			array = new double[getNumRow()*B.getNumCol()];
		}
		
		for(int i = 0; i < getSize(); i++)
			array[i] = c[i];

	return *this;
}

Matrix operator+(const Matrix& A, const Matrix& B)
{
	Matrix C = A;
	return C += B;
}

Matrix operator-(const Matrix& A, const Matrix& B)
{
	Matrix C = A;
	return C -= B;
}

Matrix operator*(const Matrix& A, const Matrix& B)
{
	Matrix C = A;
	return C *= B;
}
std::ostream& operator<<(std::ostream& out, const Matrix& C)
{
	int p = C.getNumCol();
	for(int i = 0; i < C.getSize(); i++)
	{
		if(i== p)
		{ 
			out << std::endl;
			p += C.getNumCol();
		}
		double * pC = C.getpArray();
		out << pC[i] << "\t";
	}
	out << std::endl;
	return out;
}

void Matrix::add(const double *value)
{
	for(int i = 0; i < getSize(); i++)
		array[i] = value[i];			
}

void Matrix::inverse()	//Необходимо тестирование, функция не оптимизирована
{
	assert(row == column);
	
	double E[getSize()];
	for(int i = 0; i < getSize(); i++)
		E[i] = 0;
	for(int i = 0, p = 0; i < row; i++, p += row)
		E[i+p] = 1;		
//1		
	for(int i = 0, p = 0; i < getSize(); i += row, p++)
	{
		//swap(i, max(i));
		for(int k = i + row; k < getSize(); k +=row)
		{
			double t = array[k+p]/array[i+p];
			
			for(int j = row-1; j >= 0; j--)
			{
			 	E[k+j] -=  E[i+j]*t;
				array[k+j] -=array[i+j]*t;				
			}
		}
	}
//2	
	for(int i = getSize()-1, p = 0; i > 0; i -= row, p++)
		for(int k = i - row; k >= 0; k -= row)
		{
			double t = array[k-p]/array[i-p];
			for(int j = 0; j < row; j++)
			{
				E[k-j] -= E[i-j]*t;
				array[k-p] -= array[i-p]*t;
			}
		}
//3
	for(int i = 0, p = row, j = 0; i < getSize(); i++)
	{
		if(i == p)
		{
			p += row;
			j++;
		}
		if(array[j+p-row] != 1)
			E[i] /= array[j+p-row];
	}
		
	for(int i = 0; i < getSize(); i ++)
		array[i] = E[i];	
}

void Matrix::swap(double& index1, double& index2)	//Функция не дописана
{
	double tmp = index1;
	index1 = index2;
	index2 = tmp;
	
}

int Matrix::max(const int position)
{
	if(position == getSize() -1)
		return position;
		
	double max = array[position];
	int pmax = position;
	
	for(int i = position + 1; i < getSize(); i++)
		if(fabs(max) < fabs(array[i]))
		{
			max = array[i];
			pmax = i;
		}
	
	return pmax;
}

Matrix  Matrix::transpose()
{
	double t[getSize()];
	
	for(int i=0; i < row; i++)
    	for (int j = 0; j < column; j++)
    {
        t[i*column+j] = array[j*row+i];
    }
    
	Matrix e(column,row);
	e.add(t);
//	delete []t;
	return e;
}
